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1 .  INTRODUCTION 

Many  problems  in  science  and  engineering  concern  the  determination  of 
steady-state  equilibrium  solutions  of  nonlinear  equations.  In  general,  for 
such  problems,  interest  centers  not  so  much  of  the  computation  of  a  few  specific 
equilibria  rather  than  on  an  assessment  of  the  response  of  the  system  to  the 
action  of  various  external  or  internal  influence  quantities.  In  other  words, 
we  are  interested  in  the  effect  of  changes  of  the  values  of  certain  parameters 
upon  the  computed  equilibria.  Thus,  other  than  in  the  typical  linear  case, 
for  nonlinear  problems  we  usually  have  to  consider  equations  of  the  form 

'  ••  U *  '■ 

F(z,X)  =  0  (1.1 


which  depend  nonlinearly  not  only  on  the  state  variable  z  but  also  on  a 
parameter  vector  x.  Typically,  z  varies  in  some  infinite-dimensional 

r  ■■  *  * 

space  Z  while  X  belongs  to  a  space  A  with  some  finite  dimension  m  >_1. 

In  line  with  the  indicated  objectives  it  is  not  enough  to  find  solutions 
z  of  (1.1)  for  a  few  a  priori  specified  parameter  vectors  X.  Instead,  we 
have  to  look  at  the  solutions  of  (1.1)  as  points  (z,X)  in  the  product 
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X  =  Z  x  a  of  the  state  and  parameter  space.  Under  fairly  general  conditions, 
the  set  of  all  solutions  (z,\)  of  (1.1)  in  X  forms  a  smooth  surface  --or 
more  precisely  an  m-dimensional  differentiable  manifold  --  in  that  space. 

When  (1.1)  represents  the  equilibrium  equation  of  a  mechanical  system,  this  ~ 

manifold  has  been  called  the  equilibrium  surface  of  the  system. (see  e.g.  [31]). 

t  - 

Broadly  speaking,  our  task  then  is  the  computational  analysis  of  this  -S 

solution  manifold  of  (1.1).  The  specific  aims  of  that  analysis  depend  on 
the  problem  at  hand.  For  example,  practical  needs  often  demand  the  explicit  vj 

computation  of  a  sufficient  number  of  points  in  a  selected  portion  of  the 
manifold  which  may  be  used,  for  instance,  as  the  basis  of  certain  post¬ 
processing  calculations.  Another  need  may  concern  the  determination  of  the 
stability  properties  of  the  system.  In  terms  of  the  manifold  this  leads 
typically  to  the  study  of  structural  stability  as  defined  by  modern  bifurcation 
theory.  In  essence,  an  equilibrium  point  on  the  manifold  is  structurally 
stable  if  in  some  neighborhood  all  points  on  the  manifold  have  the  same 
qualitative  behavior  and  the  task  is  to  determine  the  points  where  structural  im 

stability  is  lost  and  what  behavior  may  be  expected  there.  (See  e.g.  [33] 
for  some  interesting  applications  in  mechanical  engineering.) 

These  are  only  two  of  the  many  different  aspects  of  an  analysis  of  the 
solution  manifold  of  a  parametrized  equation  (1.1).  In  this  presentation, 
we  consider  some  of  the  questions  relating  to  a  computational  determination 
of  parts  of  the  manifold.  As  with  all  engineering  computations,  the  aim  — 

here  is  to  obtain  solutions  which  are  sufficiently  accurate  and  reliable  to 
allow  for  a  decision  about  the  system  under  study.  For  this  it  has  become 
widely  accepted  that  the  computational  procedures  should  include,  at  least,  ^ 

(i)  facilities  for  the  efficient  and  reliable  estimation  of  the  v> 

errors  of  the  computed  results;  and  -v 
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(ii)  adaptive  controls  of  the  computation  to  achieve  the  desired 
error  tolerances  with  minimal  cost. 

For  the  solution  of  certain  classes  of  linear  problems  by  finite  element 
techniques  algorithms  for  these  and  other  desirable  features  are  now 
beginning  to  be  well  understood  (see  e.g.  [1],  [10]).  But  in  the  nonlinear 
case  much  still  needs  to  be  done  before  satisfactory  procedures  of  this  type 
are  available.  The  a^m  of  this  paper  is  to  survey  some  recent  results  in  the 
area  and  to  point  to  various  open  questions. 

It  turns  out  that  a  central  aspect  of  any  computational  study  of  a  mani¬ 
fold  is  the  availability  of  simple,  but  effective  local  coordinate  systems. 
Some  approach  to  this  is  sketched  in  Section  2  below.  Then  in  Section  3  we 
turn  to  the  question  of  estimating  the  discretization  errors  between  the 
original  manifold  and  that  of  a  corresponding  discretized  equation.  It  is 
important  to  note  that  these  errors  also  depend  on  the  choice  of  the  local 
coordinate  system.  Section  4  then  presents  a  new  approach  to  the  calculation 
of  a  posteriori  estimates  of  the  discretization  errors  and  shows  their 
efficiency  in  the  case  of  a  two-dimensional  nonlinear  boundary  value  problem. 
In  Section  5  we  discuss  the  basic  form  of  a  continuation  method  for  approxi¬ 
mating  one-dimensional  manifolds  and  relate  it  to  a  general  form  of  feedback 
process.  This  in  turn  opens  up  questions  about  measuring  the  performance  of 
the  method.  Finally,  in  Section  6  we  outline  an  algorithm  for  the  adaptive 
control  of  the  discretization  errors  during  the  course  of  the  continuation 
procedure  and  show  its  performance  in  the  case  of  a  model  problem. 
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2.  LOCAL  COORDINATES 

As  noted  before,  under  certain  conditions  the  solution  manifold  of 
equation  (1.1)  is  a  finite-dimensional,  differentiable  manifold.  In  practical 
applications  (1.1)  is  usually  a  boundary  value  problem.  Thus  it  is  natural 
to  assume  that  the  mapping  F  is  a  Fredholm  mapping  of  class  Cr,  r  >_  1 , 
and  index  m  >_  1  from  an  open  subset  S  of  a  real  Banach  space  X  into 
another  such  space  Y.  We  denote  by  DF(x)  the  (Frechet)  derivative  of  F 
at  the  point  x  e  S.  A  point  x  e  S  is  a  regular  point  of  F  if  DF(x) 
maps  X  onto  Y,  and  a  point  y  e  Y  is  a  regular  val ue  of  F  if  the 
solution  set 


Fl“  (y)  =  {x  e  S,  F(x)  =  y) 


contains  only  regular  points.  Then  for  any  regular  value  y  e  Y  the  set 
/  1  \  r 

M  =  My  =  F'  '(y)  is  an  m-dimensional  C  -manifold  in  X  without  boundary 
(see  e.g.  [14]). 

For  all  computations  on  such  a  manifold  M  we  need  an  efficient  compu¬ 
tational  scheme  for  fixing  local  coordinate  systems  at  any  point  of  M. 

Among  the  many  possibilities,  a  simple  linear  scheme  appears  to  be  most 

advantageous.  In  order  to  motivate  the  idea,  let  M  be  a  one-dimensional 
13  o 

C  -manifold  in  R  .  At  a  given  point  x°  e  M  we  choose  a  one-dimensional 

3 

subspace  T  and  a  two-dimensional  plane  WcR  ,  as  shown  in  Figure  1,  such 
that  for  sufficiently  small  t  e  T  the  plane  x'  +  t  +  W  intersects  M  in 
a  unique  point  which  we  call  x(t).  In  this  way,  we  obtain  a  mapping 
t  -*•  x(t)  from  some  (small)  open  neighborhood  of  the  origin  in  R1  onto  a 
neighborhood  of  x°  in  M.  This  will  be  the  desired  local  coordinate  system 
of  M  at  x°. 


— x°+i+  & 


Figure  1 


Clearly,  T  and  W  cannot  be  chosen  arbitrarily.  Certainly  the  direct 

3 

sum  T©W  should  be  all  of  R  and,  in  order  to  get  unique  intersections 
x(t),  we  should  avoid  that  W  contains  the  tangent  manifold  ker  DF(x°)  of 
M  at  x°.  If  D  F(x°)  denotes  the  partial  derivative  of  F  at  x°  with 

W 

respect  to  W,  then  this  means  that  we  should  enforce  ker  D  F(x°)  =  {0}. 

w 

The  approach  carries  over  directly  to  the  general  case.  With  a  slight 
notational  change  the  following  result  was  proved  in  [15]: 

Theorem  2.1 :  Let  y  be  a  regular  value  of  F  and  x°  e  Suppose  that 
T,  W  are  subspaces  of  X  such  that  X  =  T  ©  W,  dim  T  =  m,  and 
ker  DwD(x°)  =  {0}.  Set  x°  =  t°  +  w°,  t°  e  T,  w°  c  T.  Then  there  exist 
open  neighborhoods  V  c  T  of  t°  and  U  c  X  of  x°  and  a  unique 
Cr-function  w:  V  -*■  W  such  that  w(t°)  =  w°  and 


a 


M  H  U  =  {x  e  X;  x=t+  w(t),  t  e  V) 


V. 


This  result  means  that  a  basis  for  the  m-dimensional  space  T  may 
serve  as  a  local  coordinate  system  on  M  near  x°.  Of  course,  if  our 
original  decomposition  of  X  into  the  state  space  Z  and  parameter  space 
A  satisfies  ker  D2F(x°)  =  {0}  then  A  itself  may  be  used  to  parametrize 
M  near  x°.  Generally,  any  x°  e  S,  where  ker  D2F(x°)  =  {0},  is  called 
a  non-singular  point  of  F  with  respect  to  the  original  (natural)  parametri 
zation;  otherwise  x°  is  a  singular  point.  In  connection  with  mechanical 
problems  the  most  frequently  occurring  singular  points  are  turning  points 
and  simple  bifurcation  points  (buckling  points).  Generally,  if  x°  e  S  is 
a  nonsingular  point,  then  the  partial  derivative  DzF(x°)  is  an  isomorphism 
of  the  state  space  Z  onto  the  range  space  Y. 


3.  ERROR  ESTIMATION 


We  consider  a  nonlinear  mapping  F:  Scx+Y  which  satisfies  the 
conditions  discussed  in  Section  2  and  assume  that  a  natural  decomposition  of 
X  =  Z  ©A  into  a  state  space  Z  and  m-dimensional  parameter  space  A  has 
been  given.  If  0  e  Y  is  a  regular  value  then  we  know  that  the  solution 
set  of  the  equation  (1.1)  is  indeed  an  m-dimensional  C  -manifold  M  in  X. 

As  noted  before,  in  applications  the  equation  (1.1)  usually  represents 
some  boundary  value  problem.  Hence  for  the  computation,  it  becomes  necessary 
to  replace  (1.1)  by  some  finite-dimensional  (discretized)  approximating 
equation.  But  since  the  parameter  space  A  is  already  finite  dimensional, 
it  is  only  the  state  variable  z  that  has  to  be  discretized.  Thus,  the 
approximating  equations,  in  general,  have  the  form 

F^(z^,x)  =  0  (3. 

where  now  F^  maps  a  discretized  space  X^  =  Z^  ©  A  to  another  such  space 

V 

If  F^  again  satisfies  the  needed  differentiability  conditions  and 
0  e  Yh  is  a  regular  value  for  it,  then  the  solution  set  of  (3.1)  forms  an 
m-dimensional  differentiable  manifold  in  X^.  Usually,  in  finite  element 
computations  we  have  X^  c.  X  and  our  task  then  is  to  assess  the  approximation 
error  between  these  manifolds  measured,  for  instance,  in  the  norm  of  X. 

The  development  of  a  rigorous  theory  of  these  errors  is  a  fairly  recent 
undertaking.  For  a  one-dimensional  parameter  space  A  a  priori  estimates 
were  first  developed  in  [13],  [17]  and  then  [14].  The  latter  results  were 
generalized  in  [15]  to  a  parameter  space  A  of  arbitrary  finite  dimension. 

All  these  results  involve  a  family  of  approximate  problems  (3.1)  which 
converge  in  some  sense  to  the  original  problem  (1.1)  when  the  real  discreti- 
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zation  index  h  >  0  tends  to  zero.  Then  it  is  shown  that  for  sufficiently 
small  h  and  in  a  neighborhood  of  some  point  of  the  solution  manifold  M  of 
(1.1)  the  approximate  problems  also  possess  solution  manifolds  for 
which  the  distance  to  M  in  X  can  be  estimated.  A  different  approach  was 
taken  in  [16].  There  only  a  single  discretized  equation  (3.1)  is  considered 
instead  of  a  converging  family  of  such  equations,  and  estimates  are  obtained 
which  are  local  in  nature  analogous  to  the  local  error  estimates  in  the 
numerical  solution  of  initial  value  problems  for  ordinary  differential  equations. 

These  a  priori  estimates  are  of  considerable  theoretical  interest.  But 
for  the  computational  task  outlined  in  the  Introduction  we  require  a  posteriori 
estimates  which  measure  the  error  of  the  specific  computed  points  on  the 
approximate  solution  manifold  M^.  However,  before  we  turn  to  these 
a  posteriori  estimates,  we  need  to  clarify  how  these  errors  are  to  be  defined. 

Without  any  further  information,  we  may  compare  a  computed  point  xh  on 
with  any  point  x°  on  M.  But  then  the  distance  ||x°-xh||x  can  hardly 
be  expected  to  represent  a  good  measure  of  the  quality  of  the  computation. 

In  order  to  correlate  points  better  we  need  to  choose  a  local  coordinate 
system  at  the  desired  exact  point  x°  on  the  manifold  M.  As  discussed  in 
the  previous  section,  this  means  that  a  coordinate  space  Tex  and 
complementary  space  W  CX  are  selected  for  which  Theorem  2.1  holds.  Then, 
when  the  approximation  is  sufficiently  good,  the  manifold  will  be  close 
to  M  and  the  same  local  coordinate  system  can  also  be  used  in  a  neighborhood 
on  Mh  of  the  point  xn  for  which  x  -xn  e  W.  This  is  schematically  shown 
in  Figure  2. 
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Thus  the  resulting  error  estimate  depends  on  the  choice  of  the  local 
coordinate  system  on  M.  This  is  an  important  point  which  cannot  be  overlooked. 
With  a  change  of  the  local  coordinate  system  the  error-norms  change.  If  x 
is  a  non-singular  point  of  M  then,  of  course,  we  may  define  the  local  co¬ 
ordinate  system  in  terms  of  the  natural  parameter  space  A  and  the  state  space 

Z.  In  other  words,  we  then  compare  x°  =  (z°,A°)  on  M  with  the  point 
h  h  o 

x  =  (z  ,A°)  on  for  which  the  parameter  vector  A0  is  the  same.  Certainly, 
at  singular  points  of  M,  such  as  the  limit  point  x  in  Figure  3,  this 
choice  fails.  Moreover,  as  the  figure  indicates,  even  at  non-singular 
points  x°  near  x  the  error-norms  ||x°-x^||^  based  on  the  natural  co¬ 
ordinate  system  may  be  unduly  large. 


4.  A  POSTERIORI  ESTIMATES 


For  the  finite-element  solution  of  various  classes  of  linear  problems 
the  theory  and  application  of  reliable,  a  posteriori  error  estimates  has 
advanced  considerably  in  recent  years  (see  e.g.  [2],  [3],  [5],  [6],  [7],  [8], 
[10]).  Not  surprisingly,  in  the  nonlinear  case  --  where  even  a  priori 
estimates  were  developed  only  recently  —  the  available  results  on  a  posteriori 
estimates  are  few.  An  interesting  approach  was  considered  in  [19],  [20],  [21], 
and  in  [10]  and  [25]  it  was  shown  that  for  certain  model  problems  in  one  space 
dimension  it  is  possible  to  generalize  the  techniques  of  the  linear  theory. 
While  the  resulting  estimates  did  prove  to  be  reliable,  their  computational 
cost  was  still  somewhat  high. 

We  sketch  here  a  new  approach  --  first  outlined  in  [27]  --  which  produces 
very  effective  and  reliable  a  posteriori  estimates  for  a  broad  class  of  non¬ 
linear  problems  and  which  turns  out  to  be  computationally  very  inexpensive. 

We  use  again  the  setting  of  the  previous  sections  and  assume  that  the 
solution  set  of  the  original  problem  (1.1)  is  an  m-dimensional ,  differentiable 
manifold  M  in  X  =  Z  x  a  and  that  represents  the  corresponding  solution 
manifold  of  some  discretized  equation  (3.1)  in  Xh  =  Zh  x  A  c  X.  As  we  saw  in 
Section  3,  any  error  estimate  depends  on  the  particular  choice  of  the  local 
coordinate  system  on  M.  The  a  posteriori  estimates  discussed  here  can  be 
developed  for  any  local  coordinate  system  for  which  Theorem  2.1  holds.  But 
in  order  to  avoid  certain  technical  details,  which  may  obscure  the  ideas,  we 
restrict  ourselves  here  to  the  case  of  a  non-singular  point  x°  =  (z°,X°)  e  M 
and  the  choice  of  the  natural  coordinate  system  based  on  the  parameter  space 
A  and  the  state  space  Z. 

We  assume  that  F  is  of  class  Cr,  r  >  2,  on  its  open  domain  S  in  X. 
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derivative  OzF(x°)  is  non-singular  and  hence  the  same  is  true  at  any  point 
x°  =  (z°,X°)  in  some  neighborhood  of  x°  in  S.  Hence  we  can  apply  a 
Newton  step  to  the  equation  G(z)  =  F(z,X°)  starting  from  z°.  The  first 
iterate  z^  is  obtained  from  the  solution  w  of  the  linearized  equation 


F(z°,X°)  +  0  F(zu,Xu)w  =  0 


=o  ,o. 


(4.1) 


by  setting  z^  =  z°  +  w.  By  the  Newton-Kantorovich  theorem  (see  e.g.  [23]) 
it  follows  that 


Ml  <  I  1 2°-Z°  I  I  <  2 1  |w|  I  , 


(4.2) 


provided  only  that  z°  is  within  the  attraction  ball  specified  by  that 
theorem  (see,  especially,  the  formulation  given  in  [24]). 

In  our  setting,  x°  =  (z°,X°)  e  is  the  computed  approximation  of 

the  desired  point  x°  =  (z°,X°).  Hence  the  required  error  ||z°-z°||  can 
be  estimated  in  terms  of  the  norm  [ | w ] |  of  the  solution  w  of  the 
linearized  equation  (4,1).  Of  course,  (4.1)  is  still  an  infinite-dimensional 
problem  which  can  only  be  solved  approximately.  But  here  we  turn  out  to  be 
exactly  in  the  setting  of  the  theory  of  linear  a  posteriori  error  estimates 
and  hence  can  apply  the  relevant  results. 

Rather  than  elaborate  on  the  general  theory  we  illustrate  the  technique 
on  the  following  one-dimensional  problem 


^  A(||)  +  B(z,s,X)  =0,  V  s  e  I  =  (0,1), 


(4.3) 


z(0)  =  z ( 1 )  =  0 


with  sufficiently  smooth  coefficient  functions  A,  B.  A  weak  formulation 
is 


N(z,X)v  =  [  [A(z')v'  -  B(z,s,X)v]ds  =  0,  V  v  e  H*(I) 

Jo  0 

and  the  corresponding  linearized  problem  (4.1)  has  the  form 


L[z,X]wv  =  N(z,X)v  + 


(DA(z‘ )w'v’ -  DuB(z,s,X)wv)ds  =  0, 


■'O 

V  v  e  hJ(I). 

For  simplicity,  suppose  that  piecewise  linear  elements  are  used  on  some 
mesh 


4:  0  =  sQ  <  s,  <  s2  <  ...  <  s„+,  ■  1.  n  -  n(i). 

In  other  words,  we  restrict  consideration  to  the  finite  dimensional  subspace 

K(A)  =  (z  c  hJ(I);  z(s)  =  J  x.^-U),  0  <  s  <  1} 

where  denotes  the  standard,  piecewise  linear  "hat  function"  with 
<J>. (s . )  =  6..,  i,j  =  l,...,n.  Then  the  exact  solution  z°  and  the  finite- 

*0  I J 

element  solution  z°  are  the  unique  functions  z°  e  H^(I)  and  z°  e  K( a) 
such  that  (4.4)  holds  for  all  v  in  H^(n)  and  K(a),  respectively. 

We  wish  to  estimate  the  norm  of  the  solution  w  of  (4.5)  at  x°;  that 
is,  of  the  problem 

L[i°,X°]wv  =0,  V  v  e  Hl(I).  (4.6) 


-  i 


•i 


Evidently,  the  zero-function  w  =  0  is  the  solution  of  (4.6)  when  v 
is  restricted  to  K(a).  Hence,  w  represents  the  error  of  the  finite- 
element  solution  of  (4.6),  and  our  problem  of  estimating  | |w| |  is  exactly 
the  problem  of  computing  an  a  posteriori  error  estimate  of  the  finite 
element  solution  w°  of  the  linearized  problem  (4.6). 

Since,  by  (4.2),  any  approximation  of  | |w| |  represents  an  a  posteriori 
error  estimate  of  the  finite  element  error  ||z°-z°||  of  the  nonlinear 
problem,  it  remains  to  apply  the  theory  of  linear  a  posteriori  estimates  to 
the  problem  (4.6).  In  the  case  of  our  model  problem,  a  simple  approach  is 
to  use  a  quadratic  finite  element  approximation  w  =  w(s),  0  s  1 ,  of 

the  solution  of  (4.6);  that  is, 


w(s)  =  piwi(s)  for  sil  <  s  <  s. 


-  /  \  (s-si-l>(si's) 

w.(s)  =  4  - — - 7j —  ,  i  =  1 ,. . .  ,n+l 


On  each  one  of  the  n+1  elements  of  A  the  evaluation  of  the  corresponding 
unknown  p..  represents  a  very  simple  calculation.  Hence,  say,  the  quantities 


(4.7) 


si 

(f  (p.-w'.(s))2ds)1/2,  i  =  1 , . . .  ,n+l 
.  11 


(4.8) 


represent  error  indicators  on  these  elements  and  our  a  posteriori  estimate 
becomes 


E  =  (  l  ni ) 
i=l  1 


2,1/2 


(4.9) 


Other  forms  of  the  error  indicators  and  other  norms  may  be  used  as  well. 

For  more  details  and  proofs  we  refer  to  [3],  [5],  [7],  [8]. 

We  shall  return  to  the  model  problem  (4.4)  in  Section  6.  As  an  example 
of  the  reliability  of  the  estimates  developed  here  we  consider  the  two- 
dimensional  problem 


N(z 


»A)v  =  (  [Aj (zs*zt)vs  +A2(2s‘zt)vt  "  C(s,t,A)v]dsdt  =  0 


V  v  E  hJ(Q)  (4.10) 


p 

where  ft  is  the  unit  square  [0,1]  x  [0,1]  in  R  .  More  specifically,  we 
use 

WZt>  =  a (zs+2t)ZS 
W*t>  *  a(zf+z^)zt 
•  to)  *  1  -  -K 

2+o 

and  choose  the  coefficient  function  C  such  that  the  exact  solution  has  the 
form 


z°(s,t)  =  A 


_ sp-s)t(l-t) _ 

Y  +  (s-0.75)2  +  (t-0.25)2 


with  some  constant  y  >  0. 


(4.11) 
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where 


M(s.t) 


D1A1  (2s*zt*  D2Al^zs’zt^ 


D1A2(zs.zt)  ^VV2^ 


k*> 


K> 


and  D.A.,  i,j  =  1,2,  denotes  the  derivative  of  the  coefficient  function 

*  J 

A.  with  respect  to  its  i-th  variable.  The  problem  (4.12)  has  the  form 
used  in  the  FEARS-system  (see  [36])  and  hence  for  the  computation  of  ||w|| 
we  can  apply  the  a  posteriori  error  estimates  developed  for  that  class  of 
linear  problems  (see  e.g.  [3],  [22]). 

The  tables  below  present  some  computational  results  obtained  in  this 
case.  Bilinear  elements  on  square,  uniform  meshes  with  the  indicated  mesh 
size  h  and  number  of  degrees  of  freedom  n  were  used  throughout.  The 
continuation  process  was  applied  to  the  one-dimensional  solution  manifold 
passing  through  the  origin  for  X  =  0.  The  points  with  parameter  values 
X  *  1,2, 3,4, 5  were  chosen  as  target  points.  In  all  cases  the  energy  norm 
for  (4.12)  was  used. 

The  table  reflects  our  general  experience  that  the  effectivity  of  the 
estimates  is  excellent.  The  performance  corresponds  to  that  reported  for 
the  linear  case  (loc.cit.) 

The  solution  (4.11)  has  a  singularity  at  (0.75,0.25)  when  the  constant 
Y  is  zero.  Thus  when  y  tends  to  zero  we  expect  the  errors  to  increase. 

This  is  indeed  the  case  but  the  performance  of  the  error  estimation  remains 
excellent  for  errors  up  to  about  10-15%.  In  Table  2  we  give  only  some  results 
for  y  =  0.1  and  the  case  h  =  1/16;  that  is,  n  =  225. 


2025 ( -1  ) 

. 2069 ( - 1  ) 

1013(-1) 

.101 7(-l ) 

6756(-l) 

.6763(-l ) 

.  4204 ( - 1 ) 

. 4396 ( - 1 ) 

.21 02 ( - 1 ) 

. 21 20(-l ) 

•  1401(-1) 

. 1 406(-l ) 

.6465 

.  6903 ( - 1 ) 

.7711 (-1) 

.  3445 (-1 ) 

. 3543(-l ) 

.  2296 ( -1 ) 

.2324(-l ) 

9553 

.1003 

.1171 

.  501 2 ( -1 ) 

.5231 (-1 ) 

.3341 (-1) 

. 3405 (- 1 ) 

.246 

.1329 

.1582 

.  6632 ( - 1  ) 

. 6968 ( - 1  ) 

.  4420 ( - 1  ) 

.4521 (-1) 
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5.  CONTINUATION  PROCESSES  ON  MANIFOLDS 


Suppose  again  that  M  denotes  the  m-dimensional  solution  manifold  of  a 
given  parametrized  equation  (1.1).  In  line  with  the  comments  in  the  Intro¬ 
duction  we  are  interested  in  computing  a  sufficient  number  of  points  on  M 
which  may  be  used,  for  instance,  as  the  basis  for  further  postprocessing. 

All  practically  useful  methods  for  this  purpose  compute  sequences  of  points 
along  prescribed  one-dimensional  submanifolds  W  of  M,  although  there  is 
certainly  room  for  some  different  approaches. 

There  are  various  ways  of  defining  such  submanifolds  N.  In  many  appli¬ 
cations  the  usual  technique  is  to  specify  a  parameter  combination  with  one 
degree  of  freedom.  For  instance,  in  a  structural  problem  the  parameter  vector 
may  represent  a  general  load  vector.  Then  the  chosen  parameter  combination 
may  determine  a  load  direction  while  the  remaining  degree  of  freedom  is  the 
load  intensity. 

The  restriction  to  the  submanifold  N  is  equivalent  with  the  construction 
of  a  modified  parametrized  equation  for  which  the  parameter  space  is  one¬ 
dimensional  (see  e.g.  [15]).  Thus  after  discretization  the  problem  reduces 
to  a  form 

G(z,X)  =  0  (5.1 

where  G:  Rn+^  -*■  Rn  is  a  given,  sufficiently  differentiable  function  of  n+1 
variables  and  n  components. 

The  one-dimensional  solution  manifold  N  of  (5.1)  may  have  several 
connected  components  and  our  objective  is  to  determine  numerically  the 
component  of  N  that  contains  a  specified  starting  point  x°  z  N. 

Most  commonly  used  processes  for  accomplishing  this 
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task  belong  to  the  class  of  predictor-corrector  continuation  methods. 

Starting  from  x°  these  methods  produce  a  sequence  x°,x\x2,...  of  points 

on  A/q.  For  any  k  ^0  the  computations  involved  in  the  step  from  x  to 
k+1 

x  correspond  essentially  to  the  construction  of  the  local  coordinate 
systems  discussed  in  Section  2. 

The  situation  is  sketched  in  Figure  4.  We  calculate  first  a  suitably 

k  k  k 

oriented  tangent  vector  p  of  at  x  .  Then  a  unit  vector  t  with 

k  T  k  k 

(p  )  t  t  0  is  chosen  to  define  the  parametrization  subspace  T  =  span(t  ) 

k 

of  a  local  coordinate  system  at  x  .  Typical  choices  are 


and 


k  k 

t  =  p  (pseudo-arclength  choice) 


k  k 

t  *  e  (local  variable  choice) 


(5.2) 


where  e\...,en+^  denote  the  unit  basis  vectors  of  Rn+1 .  The  orthogonal 
complement  of  T  is  used  as  the  space  W.  Now  a  suitable  step  along  the 

.1/  1/  b 

tangent  direction  is  chosen.  This  fixes  the  predicted  point  x  =  x  +  hkp 

-k  -  k 

and  with  it  the  (n-1 )-dimensional  plane  x  +  W  that  passes  through  x  and 

k 

is  orthogonal  to  t  .  Finally,  a  corrector  iteration  is  applied  which  --  when 
-k  -k 

started  at  x  --  produces  in  x  +  W  a  sequence  of  iterates  that  converges 

k+1 

to  the  desired  next  point  x  on  WQ.  If  the  iteration  fails  a  new  attempt 
is  made  with  a  reduced  stepsize  h^. 

In  recent  years  the  literature  on  these  continuation  processes  or 
incremental  methods,  as  they  are  also  called,  has  been  growing  at  a  rapid  pace. 
Accordingly,  there  is  no  need  to  discuss  here  any  algorithmic  details.  For 
some  surveys  of  these  processes  we  refer,  for  instance,  to  [30],  [34]  and  for 
a  description  of  a  particular  general  program  to  [28],  [29]. 


At  each  step,  these  methods  use  heuristic  techniques  to  choose  various 
quantities  for  which  only  incomplete  information  is  available,  as,  for 
instance,  the  parametrization  direction  t  and  the  steplength  hk>  Thus, 
broadly  speaking,  they  may  be  called  feedback  or  adaptive  methods.  In  [26] 
a  precise  definition  of  these  terms  was  developed  which  applies  to  many  of 
the  adaptive  procedures  in  numerical  analysis,  and  in  [4]  it  was  used  to 
prove  various  results  about  adaptive  mesh-refinement  processes  for  linear 
finite  element  computations.  There  is  certainly  a  need  to  develop  such 
results  also  for  the  continuation  methods  here  under  discussion.  But,  so 
far,  no  attempt  along  this  line  appears  to  have  been  made. 

We  shall  not  go  into  the  details  of  the  theory  in  [26]  but  identify 
only  some  of  the  main  questions.  Most  of  the  processes  under  consideration 
can  be  modeled  as  a  discrete,  state-space  system.  We  denote  the  discrete 
time  set  by  IN  =  {0,1,2,...}  and  the  state  space  by  X.  A  pair  (k,x) 


with  k  e  N,  x  e  X  is  called  an  event  and  interpreted  to  mean  that  at  time 
k  the  system  is  in  state  x.  The  system  is  assumed  to  work  on  any  problem 
p  in  a  class  P.  Its  operation  is  described  by  a  state  transition  function 
t:  WxXxP-vX  which  provides  that  when  (k,x)  is  the  current  event 
then  feedback  from  the  given  problem  p  will  lead  to  a  transition  to  the  state 
x'  =  -r(k,x,p).  Thus,  for  any  problem  p  e  P  and  starting  state  xQ  e  X,  the 
system  produces  a  trajectory  of  events 


j(p » Xg )  =  { (0 , xQ ),(1, X-|), (2, Xg ^ ^ » P )  *  ^  ^  0. 


The  entire  system  £  is  specified  by  the  three  sets  IN,  X,  and  P  and  the 
function  t.  In  view  of  the  interpretation  of  the  dependence  of  t  upon 
the  problem  p  as  a  feedback  from  that  problem,  we  call  £  a  general  feed¬ 
back  system. 

As  stated  before,  many  computational  procedures  can  be  modeled  as  such 
general  feedback  systems.  Interest  then  centers  on  measuring  their  performance. 
A  fairly  general  performance  measure  for  £  may  be  defined  as  a  mapping 
y:  n(£)  -*  M  from  the  set 


fi(£)  *  {o>(p,x  );  p  e  P.  xQ  e  X} 


of  all  possible  trajectories  of  £  into  a  given  partially  ordered  set  M 
of  performance  indices.  Recall  that  in  a  partially  ordered  set  M  an  element 
m*  e  M  is  minimal  if  there  is  no  m  e  M  such  that  m  <_ m*;  that  is,  if  for 
any  m  e  M  either  m  >  m*  or  m  and  m*  are  not  comparable.  Accordingly, 
in  [26]  a  feedback  system  £  is  called  an  adaptive  system  under  the 
performance  measure  y  if  there  exists  a  non-empty  subset  Q(y)  c  ■  x  X 
such  that  u(aj(p,xQ))  is  minimal  in  M  for  any  (p,xQ)  e  Q(y).  Th.  size 


of  the  set  Q(y)  may  be  taken  as  a  measure  of  the  robustness  of  I  (see 
e.g.  [35]). 

A  simple,  albeit  very  frequently  used  type  of  performance  measure 
distinguishes  only  between  acceptable  and  unacceptable  performance.  Here 
we  use  M  =  {0,1}  with  the  natural  ordering  and  call  the  trajectory 
co  =  co(p,xQ)  acceptable  or  unacceptable  if  p(co)  =  0  or  y(oj)  =  1, 
respectively.  Hence,  the  feedback  system  is  adaptive  under  y  if  the  set 

Q(y)  =  {(psx0)  e  P  X  x,  y(co(p,XQ) )  =  0} 
turns  out  to  be  non-empty. 

Clearly,  the  continuation  processes  discussed  above  are  feedback  systems 
under  this  definition.  But, surprisingly,  there  is  little  discussion  in  the 
literature  how  to  specify  suitable  performance  measures  for  such  processes, 
let  alone  how  to  prove  that  the  methods  are  adaptive  for  any  one  of  them. 

There  are,  of  course,  many  possible  performance  measures  that  might  be 
considered  for  continuation  methods.  If  a  specific  target  points  x*  is  to 
be  reached  on  the  given  one-dimensional  manifold  N,  then  we  may  define  y 
as  a  measure  of  the  total  number  of  points  needed  to  step  from  x°  to  x* 
or  of  the  total  work  involved  in  reaching  the  target.  But  in  many  appli¬ 
cations  no  such  target  point  is  desired  or  even  meaningful.  In  that  case, 
a  possible  choice  of  a  performance  measure  might  be  a  "velocity"  that 
indicates  how  fast  the  process  is  moving  along  Nq.  Already  simplified  model 
studies  provide  here  sometimes  rather  startling  results.  For  example, 
suppose  that  the  corrector  iteration  is  guaranteed  to  converge  whenever  the 
predicted  point  is  within  a  prescribed  distance  e  >  0  of  WQ.  Then  for 
various  model  curves  as  Nq  one  can  derive  asymptotic  relationships  for  the 
distance  s  between  successive  steps.  For  instance,  when  the  pseudo  arc- 


length  choice  (5.2)  is  used  for  the  local  parametrization,  then,  under 
certain  assumptions,  it  can  be  shown  that 


s  =  O(-)  as  c  0 

K 

where  <  is  the  curvature  of  N  .  On  the  other  hand,  when  the  local 
variable  choice  (5.3)  is  used  then,  in  general,  it  follows  that 

s  =  0(e)  as  e  ->  0. 


Figure  5 

Hence,  it  should  be  expected  that  along  strongly  curved  segments  of  the 
pseudo-arclength  process  will  take  many  small  steps  while  there  is  no  effect 
upon  the  local-variable  method.  Indeed,  this  can  be  observed  in  practical 


situations . 


These  represent  only  some  early  results  in  a  more  detailed  study  of 
suitable  performance  measures  for  continuation  methods.  But  they  indicate 
that  there  is  certainly  a  need  for  systematic  studies  in  this  area  which 
can  provide  a  basis  for  an  assessment  of  the  efficiency  of  the  various 
techniques  utilized  in  practical  programs. 


6.  ADAPTIVE  ERROR  CONTROLS 


When  a  continuation  process  is  used  for  problems  where  the  error  estimates 
of  Section  4  are  applicable,  then,  in  most  cases,  the  discretization  errors 
show  marked  variations  along  any  computed  one-dimensional  sub-manifold.  In 
fact,  it  may  even  happen  that  there  are  solutions  of  the  discretized  problem 
which  do  not  approximate  exact  solutions.  A  simple  example  of  such  spurious 
solutions  arises,  for  instance,  in  the  case  of  the  classical  Euler-Bernoull i 
rod  problem  (see  e.g.  [9]). 

As  noted  in  the  Introduction,  for  our  problems  it  is  certainly  desirable 

to  provide  facilities  in  the  continuation  procedures  that  can  control  the 

discretization  errors  along  the  computed  submanifolds.  The  basic  control 

mechanism  is  here  the  adaptive  refinement  or  de- refinement  of  the  mesh  used 

in  the  discretization.  The  continuation  process  is  applied  in  its  normal 

form  but  at  each  computed  point  a  posteriori  error  estimates  are  determined. 

k 

Then,  an  appropriate  procedure  decides  when  the  errors  at  any  point  x  are 

unacceptable.  In  that  case  a  new  mesh  A+  is  created  from  the  mesh  a  used 

in  the  computation  of  x  by  uniformly  subdividing  some  elements  of  a  and 

k  + 

coalescing  certain  contiguous  clusters  of  others.  Then  a  new  point  (x  ) 
is  determined  by  interpolation  and  the  continuation  process  for  the  new 
discretized  problem  is  started  from  (x  )  . 

Certainly,  this  type  of  combined  continuation  and  mesh-refinement 
procedure  is  a  feedback  process  in  the  sense  of  the  previous  section.  The 
question  is  how  to  define  the  decision  process  and  the  desired  performance 
measures  so  that  the  procedure  indeed  becomes  adaptive.  This  is  as  yet  an 
open  problem.  We  shall  outline  here  only  some  possible  approaches  that 
appear  to  be  very  promising. 

Our  primary  aim  is  to  compute  a  sequence  of  points  x°,x\...  which 


approximate  points  on  the  desired  manifold  nq  of  the  original  problem  in 
such  a  way  that  the  error  estimates  e°,e\...  at  these  points  satisfy 


•  ek  £  tolk  ~  6abs  +  6rel I 'x  1 1 ’  k  =  0>1’--- 


with  given  tolerance  coefficients  6  ^5  and  S^.  The  secondary  aim  is 
to  avoid  frequent  re-meshing  and  re-starting.  More  ideally  the  decision 
procedure  should  be  such  that  the  total  work  required  for  computing  a  sufficient 
number  of  points  along  a  prescribed  segment  of  S1Q  is  less  than  that  required 
for  other  choices  of  re-meshing  points  and  associated  meshes  for  which  (6.1) 
is  satisfied. 

We  restrict  this  discussion  to  the  case  of  the  problem  (4.3).  The 
corresponding  linearized  problem  (4.4)  has  the  general  form 

f  [a(s)w'v‘  -b(s)wv]ds  =  [  c(s)vds,  V  v  e  h](I)  (6. 

Jo  Jo  0 

considered  in  [7].  Any  continuous,  strictly  monotonically  increasing  function 
£  on  I  =  [0,1]  with  £(0)  =0,  £(1)  =  1,  together  with  the  specifications 


£(Sj)  ■  j  -  0,1 .  ,N=n+l 


defines  a  family  of  meshes  AU.n),  n  =  0,1,...,  on  I.  Let  w  e  H^U)  be 
the  solution  of  (6.2)  and  consider  the  partition  function 


C„(s>  -  Y„  f  [a(s)(w"(s))2]1/3ds, 


S  e  I , 


r  =  (f1  [a(s)(w"(s))2]1/3ds}_1. 


n  V 


Then  it  was  shown  in  [7]  that  under  suitable  conditions  the  family  of  meshes 
A(?0»n )  generated  by  £0  is  asymptotically  optimal  with 


t I  |e 1 1?  =  - I— 5-3  (l+0(h))  as  h  -*■  0  (6.4) 

1  E  12(n+l  )  y3 

where  h  =  max(s.-s.  j  *  1,...,N).  In  other  words,  for  any  other 

J  J  ' 

admissible  partition  function  £  the  error  for  the  mesh  A(5,n)  is  not  less 
than  that  for  A(£0,n)  for  all  sufficiently  small  h. 

Moreover,  it  was  shown  in  [7]  that  for  any  mesh  a  the  error  under  the 
energy  norm  can  be  expressed  as 


where 


o  "  1  9 

fllelll-  l  (l+0(h)) 
j-1  J 


as  h  ■*■  0 


[a(s)(w"(s))2]^3ds  (1+0(1))  as  h  -*■  0 


(6.5) 


(6.6) 


This  suggests  the  definition  of  the  piecewise  constant  function  ip  on  I 
for  which 


4»(s) 


<  s  < 


(6.7) 


Evidently,  the  value  of  ip  on  each  subinterval  represents  some  average  of 
2 

a(w")  on  that  interval,  and  we  have 


[a(s)(w“(s))2]1/3ds 


rsJ 


S  . 


j-1 


J 


^1/3(s)ds  (1+0(1))  as  h  -  0. 


(6.8) 


If  we  use  instead  of  the  e  j  in  (6.7)  the  computed  error  indicators, 
such  as  (4.8),  then  we  obtain  an  approximation  p  of  <p  and  with  it,  from 
(6.8)  and  (6.3),  an  approximate  optimal  partition  function  This  provides 

k 

us  with  a  basic  mesh- refinement  strategy.  A  point  x  is  declared  un- 

k 

acceptable  when  e  >  r  tol^  where  tol^  is  the  tolerance  defined  by  (6.1) 

1/ 

and  x  is  some  factor,  say,  t  =  0.75.  From  the  error  indicators  at  x 
the  approximate  optimal  partition  function  and  its  associated  factor 
yQ  (see  (6.3a)  and  (6.3b))  can  be  computed.  Then  (6.4)  suggests  the  use 
of  the  relation 


T  £ 


k 


12(R+1)2  Yq 


for  obtaining  an  "ideal"  mesh  size  n,  and  with  it  an  associated  mesh 
A(f^,n).  Now  the  desired  new  mesh  A+  at  xk  can  be  generated  as  an  approxi¬ 
mation  of  A(|^.n);  that  is,  we  subdivide  certain  intervals  and  re-combine 
others  so  as  to  match  A(f^\n)  as  closely  as  possible.  This  can  be 
accomplished  in  various  ways  and  we  shall  not  elaborate  upon  these  procedures. 
But  note  that  in  general  the  size  n(A+)  of  the  new  mesh  need  not  be  exactly 
equal  to  n. 

As  an  example  for  the  operation  of  this  procedure  we  consider  the  simple 
model  problem  (4.4)  with  the  coefficients 


A(o)  =  a/(l+a),  B(z,s,X)  =  A 


in  which  case  the  exact  solution  is 


30 


For  growing  X  this  solution  increases  rapidly  within  a  small  interval  near 
s  =  0.  Table  3  gives  the  results  of  the  above  procedure  when  started  with 


x  =  0  at  X  =  0.  The  tolerance  was  computed  with  5 
and  the  decision  factor  was  t  =  .75. 


,  =  .0075,  6  ,  =  .05 

abs  rel 


Figure  6  illustrates  the  changing  meshes  during  the  course  of  the 
process. 

As  the  example  indicates  the  procedure  performs  as  expected.  This  is 
our  general  experience  for  a  number  of  problems  of  varying  complexity.  The 
simplicity  of  the  approach  certainly  makes  it  very  attractive.  At  the  same 
time,  there  are  several  aspects  that  may  be  worthy  of  improvement: 

(i)  The  function  ip  used  in  the  definition  of  the  approximate, 
optimal  partition  function  £0  incorporates  only  information 
at  the  current  point  and  does  not  attempt  to  predict  the  future 
course  of  the  process. 

(ii)  The  decision  factor  t  upon  which  the  acceptance  or  rejection 
of  a  point  is  based  and  which  features  in  the  determination  of 
the  "ideal"  mesh  A(l0»n)  does  not  change  with  the  conditions 
during  the  computation. 

(iii)  In  its  present  form,  it  is  not  clear  how  to  generalize  the 
procedure  to  problems  with  higher  space  dimension. 


Preliminary  results  indicate  that  the  procedure  can  be  improved  to  account 
for  these  points  by  utilizing  some  of  the  concepts  and  approaches  developed  in 
[11]  and  [12]  for  the  case  of  parabolic  equations.  A  prototype  software 
system,  for  the  adaptive  finite  element  solution  of  a  class  of  two-dimensional, 
parametrized  non-linear  problems  is  now  under  construction  which  will 
incorporate  such  a  more  general  decision  procedure.  An  outline  of  this 
system  --  called  NFEARS  --  has  been  given  in  [32]. 
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Table  3 
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X 

^used 

Nideal 

error 

estim 

c/ 

n 

exact 

error 

% 

tol  erance 

decision 

0 

32 

0 

0 

proceed 

0.110 

32 

5 

0.1 005(-2) 

0.1361 (-2) 

0.91 02 ( -2 ) 

proceed 

0.444 

32 

13 

0.4205(-2) 

0. 1 01 3(-l ) 

0.1391  (-1) 

proceed 

1.337 

32 

22 

0 . 1 81 4  ( -1 ) 

0. 1 814  (-1 ) 

0 . 2736 ( - 1 ) 

proceed 

2.238 

32 

33 

0 . 5570 ( - 1 ) 

0 . 5586 ( - 1 ) 

0.4248(-l ) 

refine 

33 

24 

0.3077(-l ) 

0. 3077(-l ) 

0.4261  (-1  ) 

proceed 

2.350 

33 

26 

0 . 3456 ( - 1 ) 

0 . 3465 ( -1  ) 

0 . 4466 ( -1  ) 

proceed 

2.685 

33 

32 

0 . 4947 ( - 1  ) 

0 . 5242 ( -1  ) 

0.5105(-1) 

proceed 

3.589 

33 

39 

0.1451 

0.1466 

0.7130 

refine 

37 

27 

0.4777 (-1 ) 

Q.4776H) 

0. 71 93 ( -1 ) 

proceed 

3.700 

37 

28 

0 . 5248 ( -1 ) 

0.5251 (-1) 

0 . 7484 ( - 1 ) 

proceed 

4.036 

37 

32 

0.71 63(-l ) 

0 . 7606 ( - 1  ) 

0. 8399( -1 ) 

proceed 

4.938 

37 

45 

0.2082 

0.2094 

0.1147 

ref i ne 

45 

44 

0.1112 

0.1109 

0.1153 

proceed 

5.049 

45 

48 

0.1267 

0.1267 

0.1198 

ref i ne 

44 

31 

0. 781 0(-l ) 

0.7857(-l ) 

0.1202 

proceed 

5.161 

44 

31 

0.7620(-l) 

0. 8644 ( -1  ) 

0.1249 

proceed 
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